data {
  int<lower=1> N;                   // number of observations 
  int<lower=1> C;                   // number of congresses (time periods)
  int<lower=1,upper=C> cong[N];     // map between observations & congresses
  int<lower=1,upper=56> nVotes[N];  // total votes
  int<lower=0,upper=55> nWins[N];   // majority party victories
  real lvRatio[N];                  // log(ratio of avg vote-shares)
  matrix[C,C] Pinverse;             // inverse of penalty matrix (for GMRF prior)
}
transformed data {
  matrix[C,C] chol_Pinverse = cholesky_decompose(Pinverse);  
}
parameters {
  real<lower=0> phi_raw;    
  vector[C] lambda_raw;      
  vector[C] rho_raw;
  real<lower=0> tau_sq_lambda;  
  real<lower=0> tau_sq_rho;
}
transformed parameters {
  real<lower=0> phi;
  vector[C] lambda;      
  vector[C] rho;
  
  // implies phi ~ normal(0, 100)
  phi = 0 + 100 * phi_raw;
  
  // non-centered parameterization
  // implies lambda and rho are multivariate normal
  lambda  = (tau_sq_lambda * chol_Pinverse) * lambda_raw;
  rho  = (tau_sq_rho * chol_Pinverse) * rho_raw;
}
model {
  // local/temporary variables
  real logLik;      // log likelihood
  real logPrior;    // log prior
  vector[N] alpha;  // shape1 for Beta-Binomial 
  vector[N] beta;   // shape2 for Beta-Binomial

  // log likelihood
  for (n in 1:N) {
    real eta_n;
    real theta_n;
    eta_n = lambda[cong[n]] + rho[cong[n]] * lvRatio[n];
    theta_n   = inv_logit(eta_n);    
    alpha[n] = theta_n * phi;
    beta[n]  = (1 - theta_n) * phi;
  } 
  logLik = beta_binomial_lpmf(nWins | nVotes, alpha, beta);
  
  // log priors
  logPrior = 
    normal_lpdf(tau_sq_lambda | 0, 1) +
    normal_lpdf(tau_sq_rho | 0, 1) +
    normal_lpdf(phi_raw | 0, 1) +
    normal_lpdf(lambda_raw | 0, 1) + 
    normal_lpdf(rho_raw | 0, 1);
  
  // log posterior (up to proportion)
  target += logPrior + logLik; 
}
generated quantities {
  vector[C] bias = inv_logit(lambda) - 0.5;  // estimated bias toward majority party
}
